
#ifndef __RBFFUNCTION_3__
#define __RBFFUNCTION_3__

#include "../typedefs.h"
#include <vector>

class RbfFunction {
   
public:
   
   RbfFunction(): pPoints(0), npolyterms(0) { }
   
   virtual ~RbfFunction() { }
   
   /// It computes the RBF inverse matrix
   void compute_iM(vector<Point3> &c) {
      
      centers = c;
      const uint &ncenters = centers.size();
      
      /// Built the RBF matrix 
      unsigned int size = ncenters+npolyterms;
      Matrix M(size, size);
      for (unsigned int row=0; row<ncenters; ++row) {
         const Point3 &Pf = centers[row];
         for (unsigned int col=0; col<ncenters; ++col) {
            const Point3 &Pc = centers[col];
            M(row,col) = phi(Pf-Pc);
         }
         for (unsigned int col=ncenters; col<size; ++col) {
            M(row,col) = polynomial(Pf,col-ncenters);
         }
      }
      
      /// Fill the bottom rows of the matrix
      for (unsigned int row=ncenters; row<size; ++row) {
         for (unsigned int col=0; col<size; ++col) {
            if (col < ncenters) M(row,col) = M(col,row);
            else                M(row,col) = 0.0;
         }
      }
      
      /// Invert the matrix
      iM = i(M);
   }
   
   void setEvalPoints(vector<Point3> &points) {
      pPoints = &points;
   }
   
   Number phi(const Vector3& v) {
      const Number &sqrLen = v*v;
      //~ return sqrLen*(sqrt(sqrLen)); // phi(v)=|v|^3
      return std::sqrt(sqrLen);   // phi(v)=|v|
   }
   
   Number polynomial(const Point3& p, unsigned int i) {
      
      assert(i<10);
      Number value=0.0;
      switch(i) {
         case 0:  value=1.0;             break;
         case 1:  value=p[0];            break;
         case 2:  value=p[1];            break;
         case 3:  value=p[2];            break;
         case 4:  value=p[0]*p[1];       break;
         case 5:  value=p[0]*p[2];       break;
         case 6:  value=p[1]*p[2];       break;
         case 7:  value=p[0]*p[0];       break;
         case 8:  value=p[1]*p[1];       break;
         case 9:  value=p[2]*p[2];       break;
      }
      return value;
   }
   
protected:
   
   Matrix            iM;              // the rbf inverse matrix
   vector<Point3>    centers;         // the rbf centers
   vector<Point3>   *pPoints;         // a pointer to vector of points that will be evaluated
   uint              npolyterms;      // number of polynomial terms
};

#endif
